home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / pascal / tpl60n19.zip / ARISOURC.ZIP / F48FSIN.ASM < prev    next >
Assembly Source File  |  1993-01-24  |  8KB  |  164 lines
  1.  
  2. ; *******************************************************
  3. ; *                                                     *
  4. ; *     Turbo Pascal Runtime Library Version 6.0        *
  5. ; *     Real Sine/Cosine                                *
  6. ; *                                                     *
  7. ; *     Copyright (C) 1989-1992 Norbert Juffa           *
  8. ; *                                                     *
  9. ; *******************************************************
  10.  
  11.              TITLE   F48FSIN
  12.  
  13.              INCLUDE SE.ASM
  14.  
  15.  
  16. CODE         SEGMENT BYTE PUBLIC
  17.  
  18.              ASSUME  CS:CODE
  19.  
  20. ; Externals
  21.  
  22.              EXTRN   RealAdd:NEAR,RealSub:NEAR,RealPoly:NEAR,RealMulNoChk:NEAR,
  23.              EXTRN   HaltError:NEAR,RealTrunc:NEAR,RealMul:NEAR,RealFloat:NEAR
  24.              EXTRN   RealReduceMP:NEAR,RFrac:FAR,RealMulFNoChk:NEAR
  25.  
  26. ; Publics
  27.  
  28.              PUBLIC  RSin,RCos
  29.  
  30. ;-------------------------------------------------------------------------------
  31. ; RSin and RCos are different entries into the same routine that computes the
  32. ; sine as well as the cosine, as cosine is the same as the sine with an offset
  33. ; of pi/2. The argument is reduced to the range -pi/4..+pi/4. Depending on the
  34. ; octant into which the argument falls and the function to be computed, the
  35. ; cosine or sine of the reduced argument is computed by polynomial approxi-
  36. ; mations to the sine/cosine. To prevent excessive loss of precision in the
  37. ; argument reduction, the mapping of the argument to the reduced argument is
  38. ; done via a pseudo extended precision calculation. However this extended
  39. ; precision calculation is only available for arguments inside the range
  40. ; -1.68e+9..1.68e+9. Outside this range, a simple argument reduction is used.
  41. ; This causes a sudden drop in accuracy if the absolute value of the argument
  42. ; gets larger than 1.68e+9. For arguments bigger than 6.9e12, there is a total
  43. ; loss of precision. Outside this range, Sin always returns 0, Cos returns 1.
  44. ; The approximating polynomial for the sine on -pi/4..pi/4 is:
  45. ;
  46. ; P(x) := ((((2.476053850842e-8 *x² + 2.755533965768e-6)*x² -
  47. ;             1.984126372865e-4)*x² + 8.333333325083e-3)*x² -
  48. ;             1.666666666663e-1)*x²*x + x
  49. ;
  50. ; This approximation has a theoretical maximum relative error of 1.13e-14.
  51. ; Maximum observed error when evaluated in REAL arithmetic is 1.04e-13.
  52. ;
  53. ; The approximating polynomial for the cosine on -pi/4..pi/4 is:
  54. ;
  55. ; P(x) := ((((-2.715314622037e-7 *x² + 2.479862035332e-5)*x² -
  56. ;              1.388887879411e-3)*x² + 4.166666651281e-2)*x² -
  57. ;              4.999999999923e-1)*x² + 1.000000000000e+0
  58. ;
  59. ; This approximation has a theoretical maximum relative error of 1.41e-13.
  60. ; Maximum observed error when evaluated in REAL arithmetic is 1.49e-12.
  61. ;
  62. ; INPUT:     DX:BX:AX  argument
  63. ;
  64. ; OUTPUT:    DX:BX:AX  sin of argument
  65. ;
  66. ; DESTROYS:  AX,BX,CX,DX,SI,DI,Flags
  67. ;-------------------------------------------------------------------------------
  68.  
  69. RCos         PROC    FAR
  70.              PUSH    BP                ; save caller's frame pointer
  71.              AND     DH, 7Fh           ; abs(x), since cos (x) = cos (abs(x))
  72.              MOV     BP, 2             ; load function code for cosine
  73.              JMPS    $cos_entry        ; goto sine/cosine computation
  74. RCos         ENDP
  75.  
  76.              ALIGN   4
  77.  
  78. RSin         PROC    FAR
  79.              PUSH    BP                ; save caller's frame pointer
  80.              XOR     BP, BP            ; load function code for sine
  81. $cos_entry:  MOV     CH, 80h           ; load mask to extract sign bit
  82.              AND     CH, DH            ; extract sign bit
  83.              XOR     DH, CH            ; x := abs(x)
  84.              PUSH    CX                ; save sign
  85. $value_ok:   MOV     CX, 04E81h        ; load
  86.              MOV     SI, 0836Eh        ;  real constant
  87.              MOV     DI, 022F9h        ;   4/pi
  88.              CMP     AL, 9Fh           ; x >= 2^30 ?
  89.              JAE     $big_val          ; x/(pi/4) too big to be stored in long
  90.              PUSH    AX                ; save x
  91.              PUSH    BX                ;  on
  92.              PUSH    DX                ;   stack
  93.              CALL    RealMulFNoChk     ; compute x/(pi/4)
  94.              XOR     CH, CH            ; signal truncation desired
  95.              CALL    RealTrunc         ; n = Trunc (x/(pi/4))
  96.              ROR     AL, 1             ; if quadrant odd, set CF=1
  97.              ROL     AL, 1             ; restore register, CF=1 if odd quadrant
  98.              ADC     AX, 0             ; increment quadrant
  99.              ADC     DX, 0             ;  if quadrant odd
  100.              ADD     BP, AX            ; q := quadrant + flag
  101.              CALL    RealFloat         ; convert quadrant to REAL value
  102.              POP     DI                ; get
  103.              POP     SI                ;  back
  104.              POP     CX                ;   x
  105.              TEST    BP, 3             ; determine octants that use cos,sin
  106.              PUSHF                     ; save sin/cos flag
  107.              PUSH    BP                ; save q
  108.              MOV     BP, OFFSET c1     ; pointer to reduction constant table
  109.              CALL    RealReduceMP      ; compute x + n*(pi/2) to full precision
  110.              MOV     CX, 5             ; use five coefficients
  111.              XOR     SI, SI            ; signal P(x²)*x+x wanted
  112.              POP     BP                ; get back q
  113.              POPF                      ; get back sin/cos flag
  114.              JPE     $sine             ; octants 0, 3, 4, 7 take sine
  115. $cosine:     MOV     DI,OFFSET CS_COEFF; pointer to 1st coeff. of cos approx.
  116.              INC     SI                ; signal P(x²) wanted
  117.              CALL    RealPoly          ; P(x²) in DX:BX:AX
  118.              MOV     CX, 00081h        ; load
  119.              XOR     SI, SI            ;  floating-point
  120.              MOV     DI, SI            ;   constant 1.0
  121.              CALL    RealAdd           ; 1 + P(x²)
  122.              JMPS    $make_sign        ; put appropriate sign on result
  123. $sine:       MOV     DI,OFFSET SN_COEFF; pointer to 1st coeff. of sin approx.
  124.              CALL    RealPoly          ; P(x²)*x+x in DX:BX:AX
  125. $make_sign:  POP     CX                ; get sign mask
  126.              INC     BP                ; determine if
  127.              TEST    BP, 4             ;  result has to negated
  128.              JZ      $ende             ; result ok
  129.              XOR     CH, 80h           ; negate sign flag
  130. $ende:       POP     BP                ; restore caller's frame pointer
  131.              OR      AL, AL            ; result zero ?
  132.              JZ      $end_sin          ; yes, done
  133.              XOR     DH, CH            ; make sign bit
  134. $end_sin:    RET                       ; done
  135.  
  136. $big_val:    SUB     CL, 3             ; 1/(2*pi) in DI:SI:CX
  137.              CALL    RealMulNoChk      ; x/(2*pi)
  138.              CALL    RFrac             ; frac (x/(2*pi))
  139.              MOV     CX, 02183h        ; load
  140.              MOV     SI, 0DAA2h        ;  constant
  141.              MOV     DI, 0490Fh        ;   2*pi
  142.              CALL    RealMulNoChk      ; 2*pi * frac (x/(2*pi)) = remainder
  143.              JMP     $value_ok         ; reduction to 0..2*pi done
  144.  
  145. SN_COEFF     DB      067h,               0B1h,0D4h  ; -2.476053850842e-8
  146.              DB      06Eh,05Ch,08Bh,0B6h,0EBh,038h  ;  2.755533965768e-6
  147.              DB      074h,0B5h,09Ch,0FCh,00Ch,0D0h  ; -1.984126372865e-4
  148.              DB      07Ah,044h,086h,088h,088h,008h  ;  8.333333325083e-3
  149.              DB      07Eh,0A9h,0AAh,0AAh,0AAh,0AAh  ; -1.666666666663e-1
  150. CS_COEFF     DB      06Bh,               0C7h,091h  ; -2.715314622037e-7
  151.              DB      071h,034h,0B7h,0A1h,006h,050h  ;  2.479862035332e-5
  152.              DB      077h,02Eh,00Ah,058h,00Bh,0B6h  ; -1.388887879411e-3
  153.              DB      07Ch,018h,0A0h,0AAh,0AAh,02Ah  ;  4.166666651281e-2
  154.              DB      07Fh,0EFh,0FFh,0FFh,0FFh,0FFh  ; -4.999999999923e-1
  155.  
  156. C1           DB      080h,000h,000h,000h,00Fh,0C9h  ; pi/4-eps =
  157. C2           DB      070h,0c2h,068h,021h,0a2h,0dah  ;      eps =
  158.  
  159.              ALIGN   4
  160.  
  161. CODE         ENDS
  162.